No description has been provided for this image

Do We Really Need to Remove Euler Deconvolution Solutions Based on Spatial or Other Criteria Besides Residuals in Magnetic and Gravity Data?¶

In magnetic and gravity data interpretation, Euler deconvolution is widely used for estimating source locations and depths based on potential field data. It is a popular method due to its simplicity and ability to work with large datasets. However, like any geophysical technique, Euler deconvolution has its limitations, often yielding solutions with varying degrees of accuracy. One key question is whether we should filter these solutions based on spatial or other criteria, beyond simply relying on the residuals.

Understanding Euler Deconvolution Solutions¶

Euler deconvolution relies on the assumption of a structural index that describes the geometry of the causative body (e.g., sphere, cylinder, dyke). This technique solves for the source location and depth by minimizing the difference between the observed and modeled field using Euler’s homogeneity equation. The result is a set of solutions that include position, depth, and residuals.

Euler Deconvolution¶

Euler Deconvolution is a widely used technique in geophysics for interpreting potential field data, such as magnetic and gravity fields. It is employed to estimate the location and depth of geological bodies, such as dykes, using the observed anomalies in the potential field.

Euler's Equation¶

The fundamental concept of Euler Deconvolution is based on Euler's equation for potential fields. For a point source, the relationship between the observed potential field ( T ), the coordinates ( (x, y, z) ), and the depth ( d ) of the source is given by:

$$ T = \frac{A}{\left((x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2\right)^{\frac{n}{2}}} $$

where:

  • ( T ) is the total magnetic field or gravity anomaly.
  • $( (x_0, y_0, z_0) )$ are the coordinates of the source.
  • ( n ) is the structural index which indicates the type of geological body.
  • ( A ) is a scaling factor related to the strength of the anomaly.

Euler Deconvolution Formula¶

Euler Deconvolution uses a local approach to estimate the depth and location of the source body by solving a system of linear equations. The Euler equation is:

$$ \frac{\partial T}{\partial x} \cdot (x - x_0) + \frac{\partial T}{\partial y} \cdot (y - y_0) + \frac{\partial T}{\partial z} \cdot (z - z_0) = (T - T_0) \cdot * SI $$

where:

  • $( \frac{\partial T}{\partial x} )$, $( \frac{\partial T}{\partial y} )$, and $( \frac{\partial T}{\partial z} )$ are the gradients of the potential field in the x, y, and z directions respectively.
  • $( T_0 )$ is the background field value.

Solution Process¶

  1. Gradient Calculation: Compute the gradients $( \frac{\partial T}{\partial x} )$, $( \frac{\partial T}{\partial y} )$, and $( \frac{\partial T}{\partial z} )$ from the potential field data.

  2. System of Equations: For each grid point in the data, form the linear system of equations based on Euler’s equation. Solve these equations to obtain the estimated coordinates $( (x_0, y_0, z_0) )$ and the residuals.

  3. Residual Analysis: Analyze the residuals to evaluate the fit of the model. Lower residuals indicate a better fit.

Magnetic Dyke Depth Prediction¶

Magnetic dykes are linear geological bodies that cause anomalies in the magnetic field. Euler Deconvolution can be applied to predict the depth of these dykes by analyzing the magnetic anomalies they produce.

Depth Estimation¶

For a dyke, which is modeled as a vertical or near-vertical body, the structural index ( n ) is typically set to ( 1 ). The depth of the dyke can be estimated from the Euler Deconvolution as follows:

$$ d = \frac{A \cdot \left((x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2\right)^{\frac{1}{2}}}{T - T_0} $$

where:

  • ( d ) is the depth of the dyke.
  • ( (x, y, z) ) are the coordinates of the observation point.
  • ( T ) is the observed magnetic field anomaly.

Example Application¶

  1. Data Preparation: Collect magnetic anomaly data and compute the gradients.

  2. Deconvolution: Apply Euler Deconvolution to estimate the depth and location of the dyke.

  3. Visualization: Visualize the results using contour plots and scatter plots to assess the depth and location of the dyke.

Euler Deconvolution is a powerful tool for interpreting potential field data and estimating the location and depth of geological bodies such as magnetic dykes. By analyzing the observed anomalies and solving the associated equations, geophysicists can gain valuable insights into subsurface structures.

Challenges in Euler Deconvolution Interpretation¶

  1. Non-unique Solutions: One major issue with Euler deconvolution is that it can produce a large number of solutions, some of which may be physically unrealistic or influenced by noise.

  2. Spatial Clustering: The solutions may also cluster around areas where data noise is high, leading to spurious or misleading interpretations.

  3. Uncertainty in Structural Index: Incorrect assumptions about the structural index can lead to erroneous depth and location estimates.

  4. Depth Ambiguities: Shallow sources may be misinterpreted as deep, or vice versa, depending on the data quality and spatial sampling.

Residuals vs. Spatial or Other Criteria for Filtering Solutions¶

The residual is often used to filter Euler deconvolution solutions, under the assumption that lower residuals indicate better fits and therefore more reliable solutions. However, residuals alone may not always provide sufficient filtering. In certain cases, relying only on residuals can lead to the retention of incorrect solutions, especially in areas of high noise or complex geology. Here’s why incorporating spatial or other criteria might be necessary:

1. Spatial Continuity¶

In many geophysical scenarios, geological structures exhibit spatial continuity. Therefore, solutions that are spatially isolated from others, especially in regions where the geology suggests continuity, could be discarded. This is particularly useful in complex terrains where noise may generate isolated solutions with low residuals but no geological relevance.

2. Depth Consistency¶

Depth estimates should be consistent with the broader geologic context. If a set of solutions yields depths that are inconsistent with known geological structures (e.g., excessively shallow or deep solutions in a homogenous area), they could be removed.

3. Structural Geometry¶

If the Euler deconvolution solution yields source geometries that do not align with known or expected structural trends (e.g., lineaments, faults, or dykes), it may be worth rejecting such solutions. For example, in regions with predominantly planar structures, solutions that suggest spherical sources may be incorrect.

4. Correlation with Other Data¶

In areas where multiple geophysical datasets are available (e.g., gravity, magnetic, seismic), solutions can be filtered based on how well they correlate with other datasets. If the Euler deconvolution solutions from magnetic data do not align with gravity data interpretations, this could signal the need for rejection or further investigation.

5. Clustering and Grouping¶

Solutions that cluster within a certain depth range or spatial extent might be more reliable, whereas outliers could be rejected. Grouping methods, like k-means clustering or other statistical tools, could help in isolating meaningful groups of solutions.

Caution: Potential Loss of Information During Filtering¶

While filtering based on spatial, geological, or other criteria improves the accuracy of Euler deconvolution results, there is also a risk of losing valuable information. Some solutions that may initially appear as outliers could correspond to real, yet subtle, geological features that are not well represented in the broader dataset. Over-filtering can lead to the exclusion of these potentially significant features, especially in complex geophysical environments. Therefore, careful consideration is required to strike a balance between removing noise and preserving useful data.

Conclusion: Filtering Beyond Residuals is Important¶

While residual filtering is a necessary step in cleaning up Euler deconvolution solutions, it is not sufficient on its own. The incorporation of spatial, geological, and data-driven criteria can significantly improve the reliability of the interpretations. Filtering solutions based on spatial clustering, depth consistency, correlation with known geological structures, and coherence with other geophysical datasets provides a more holistic approach to data interpretation.

Thus, for magnetic and gravity data, the removal of Euler deconvolution solutions based on additional criteria beyond residuals can greatly enhance the accuracy and usefulness of the results, reducing the risk of misinterpretation in complex geophysical environments.

Demonstration of Euler Deconvolution on a Synthetic Model and Real Aeromagnetic Dataset¶

Euler deconvolution is a widely used technique in geophysics for estimating the depth and location of magnetic or gravity sources. To fully understand its capabilities and limitations, we will first demonstrate Euler deconvolution using a synthetic magnetic model, followed by an application to a real aeromagnetic dataset. This two-step approach allows us to evaluate the method in both controlled and real-world scenarios.

Step 1: Applying Euler Deconvolution on a Synthetic Model¶

We start by building a synthetic magnetic model designed to simulate the magnetic response from simple geological structures, such as a buried sphere or cylinder. This model is used to generate magnetic anomaly data, which mimics the behavior of real magnetic data but with known source locations and depths.

In this synthetic model, we assume several magnetic sources at specific depths beneath the surface. The magnetic anomaly generated by these sources is plotted and forms the basis of our analysis. By applying Euler deconvolution to this synthetic dataset, we can evaluate how well the technique recovers the true source positions.

Evaluation of Results on the Synthetic Model¶

The Euler deconvolution solutions are expected to match the known source depths and locations of the synthetic model. However, various factors, such as noise and incorrect assumptions about the structural index, can lead to deviations in the results. For instance, some solutions may scatter around the true source location due to data noise, while others may yield inaccurate depths if the structural index is not properly chosen.

In this controlled environment, we observe that:

  • Euler deconvolution successfully recovers the approximate depth and location of the sources when the structural index is correctly defined.
  • The presence of noise slightly affects the solutions, but spatial clustering of solutions helps in identifying the true sources.
  • Some outliers can be discarded based on spatial filtering techniques or by comparing the results with residual errors.

This synthetic test demonstrates the potential of Euler deconvolution, but also highlights the importance of correctly selecting the structural index and applying additional filtering techniques to improve the reliability of the results.

Step 2: Application of Euler Deconvolution on a Real Aeromagnetic Dataset¶

Next, we apply Euler deconvolution to a real-world aeromagnetic dataset from a geologically complex area. Aeromagnetic surveys are widely used in mineral exploration, and the dataset used here reflects magnetic anomalies caused by subsurface structures such as faults, dykes, or mineralized bodies.

The real dataset contains magnetic anomalies across a large area, and our goal is to identify potential source locations and depths. The Euler deconvolution method is applied to the gridded magnetic data, and solutions are generated for each data point.

Evaluation of Real Dataset Results¶

The results from applying Euler deconvolution to the real aeromagnetic dataset are more challenging to interpret compared to the synthetic model. Unlike the controlled synthetic data, the real dataset contains:

  • Noise from geological complexity, survey errors, and natural magnetic variations.
  • Overlapping anomalies caused by multiple geological bodies, leading to more complicated magnetic signals.
  • Ambiguities in depth estimation, as the exact structural indices for the underlying sources are not known with certainty.

To filter the solutions, we use additional criteria beyond residual minimization, such as spatial clustering and consistency with known geological features. Solutions that are spatially isolated or inconsistent with the geological context are removed, leaving us with a more reliable set of source depth estimates.

Key Observations¶

  • The application of Euler deconvolution to real aeromagnetic data produces a large number of solutions, not all of which are geologically meaningful. Filtering based on spatial and depth consistency helps in improving the reliability of the results.
  • In regions with known geological structures, the Euler deconvolution results correlate well with existing geological maps, confirming the method’s effectiveness in identifying subsurface sources.
  • Some depth estimates show discrepancies, highlighting the limitations of the method in areas with complex geology and multiple overlapping anomalies.

Through this exercise, we have demonstrated the application of Euler deconvolution on both synthetic and real-world magnetic datasets. The technique proves to be an effective tool for depth and location estimation, but care must be taken in interpreting the results, especially when applied to real-world data. Using additional filtering criteria beyond residuals—such as spatial consistency, geological context, and correlation with other geophysical datasets—greatly enhances the reliability of Euler deconvolution solutions.

In [1]:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import Forward_Model as fm
import time
from tqdm import tqdm
import grav_mag_inv as gmi
import Grav_Mag_Transform as gt
import GravMagPro as gmp
import rasterio
import Euler_Deconvolution as ev
In [2]:
length = 800
width = 800
height = 100
center_x = 500
center_y = 500
depth_z = 750
susceptibility_cuboid = 0.07
susceptibility_background = 0.0004
susceptibility1, X, Y, Z = fm.create_3d_cuboid_oriented(length, width, height, center_x, center_y, depth_z, 
                                                    susceptibility_cuboid, susceptibility_background, azimuth_angle=45, tilt_angle=0, L=1000, N=20)
# Parameters
L = 1000  # Size of the 3D grid (meters)
N = 20   # Number of points in each dimension
cylinder_radius = 50  # Radius of the cylinder (meters)
cylinder_length = 600  # Length of the cylinder along the Z-axis (meters)
cylinder_center = np.array([L/2, L/2])  # Center of the cylinder in the X-Y plane
susceptibility_cylinder = 0.05  # SI
susceptibility_background = 0.0004  # SI
susceptibility2, X, Y, Z = fm.create_3d_cylinder(cylinder_radius, cylinder_length, susceptibility_cylinder, susceptibility_background, L = L, N = N)
susceptibility3 = susceptibility1 + susceptibility2
fig = gmi.plot_3d_subsurface(susceptibility3, X, Y, Z)
fig.show()
In [3]:
# Earth's magnetic field parameters
I = np.radians(90)  # Magnetic inclination in radians
D = np.radians(0)  # Magnetic declination in radians
B_0 = 31000  # Earth's magnetic field strength in nT (assumed)
In [4]:
start_time = time.time()
delta_T_sim = fm.create_3d_forward_model_Fast(I, D, B_0, susceptibility3, L, N)
elapsed_time = time.time() - start_time
print(f"Time taken: {elapsed_time:.2f} seconds")

Time taken: 0.45 seconds
In [5]:
fm.plot_forward_model(delta_T_sim, X, Y)
No description has been provided for this image
In [6]:
T_sim = delta_T_sim - np.mean(delta_T_sim)
In [7]:
cell_size_x = L/N
cell_size_y = L/N
In [8]:
# Automatically determine the upper-left corner
upper_left_x = 1000
upper_left_y = 1000
In [9]:
# Create an affine transform for the grid (Y cell size is negative to go downward)
transform = rasterio.transform.from_origin(upper_left_x, upper_left_y, cell_size_x, -cell_size_y)
In [10]:
# Optional parameters for the Grid class
nodata_value = -9999
name = 'OptimizedGrid'
filename = 'theoritical_grid.tif'
mask = None  # No mask for now
crs = None  # Set your CRS here
In [11]:
data_array = T_sim
dMx, dMy, dMz = ev.compute_derivatives(data_array, transform, nodata_value, name, filename, mask, crs = None)
In [12]:
x_grid, y_grid = np.meshgrid(
    np.linspace(0, 1000, 20),
    np.linspace(0, 1000, 20)
)
In [13]:
ev.plot_contours(x_grid, y_grid, dMx, label = 'nT/m', title = 'Horizontal Derivative in X-direction')
No description has been provided for this image
In [14]:
ev.plot_contours(x_grid, y_grid, dMy, label = 'nT/m', title = 'Horizontal Derivative in Y-direction')
No description has been provided for this image
In [15]:
ev.plot_contours(x_grid, y_grid, dMz, label = 'nT/m', title = 'Vertical Derivative')
No description has been provided for this image
In [16]:
max_depth = 1000
window_size=3
slide_distance=1
structural_index=1 
min_denominator=1e-8
In [17]:
distance_threshold = 700
residual_threshold = 2000
# Apply Euler Deconvolution using a sliding window in X and Y directions
structural_index = 1  # Example structural index
start_time = time.time()

solutions = ev.euler_deconvolution_3d_window_xy_n(T_sim , dMx, dMy, dMz, structural_index, window_size,
                                               distance_threshold, residual_threshold)

Processing Euler Deconvolution: 100%|██████████████████████████████████████████| 324/324 [00:00<00:00, 4093.66window/s]
In [18]:
elapsed_time = time.time() - start_time
print(f"Time taken: {elapsed_time:.2f} seconds")

Time taken: 0.12 seconds
In [20]:
# Convert results to a DataFrame for easier visualization and further processing
solutions_df = pd.DataFrame(solutions)
In [21]:
solutions_df.dropna(inplace=True)
solutions_df
Out[21]:
x0 y0 z0 residuals
0 2637.911213 3441.673654 6052.772698 85.205662
1 168.359932 351.354019 424.154735 1076.856072
2 560.013633 -41.759949 -745.806500 824.035410
3 303.230882 -110.203398 -786.714813 439.542265
4 7.723873 -199.951616 -1102.907099 82.299506
... ... ... ... ...
319 694.496594 801.531202 9.032875 2.836657
320 749.145076 792.346964 -44.382645 2.124641
321 777.405247 813.386558 86.541204 0.270486
322 836.325891 811.418009 62.514197 2.935030
323 917.977263 782.943154 -149.233340 0.132659

324 rows × 4 columns

In [22]:
solutions_df = solutions_df[solutions_df['z0'] >=5]
solutions_df = solutions_df[solutions_df['x0'] >=0]
solutions_df = solutions_df[solutions_df['y0'] >=0]
solutions_df = solutions_df[solutions_df['residuals'] <=30]
In [23]:
import seaborn as sns
In [24]:
sns.histplot(data=solutions_df, x='residuals', bins=100)
plt.xlabel('Residuals')
plt.ylabel('Count')
plt.title('Histogram of Residuals')
plt.show()
No description has been provided for this image
In [25]:
sns.histplot(data=solutions_df, x='z0', bins=100, kde=True)
plt.xlabel('Residuals')
plt.ylabel('Count')
plt.title('Histogram of Depth')
plt.show()
No description has been provided for this image
In [26]:
ev.plot_results(solutions_df, T_sim, x_grid, y_grid)
No description has been provided for this image
In [27]:
df = pd.read_csv('mag_data11.csv')
df.head()
Out[27]:
X Y Z mag
0 330000 6629200 211.055369 5535
1 330050 6629200 210.654999 5500
2 330100 6629200 210.286555 5476
3 330150 6629200 209.952679 5479
4 330200 6629200 209.655742 5496
In [28]:
mag_grid, x_grid, y_grid = ev.create_grid(df, 200, 200)
In [29]:
mag_grid.shape
target_points_x = 201
target_points_y = 201
In [30]:
# Extracting data from DataFrame
x = df['X'].values
y = df['Y'].values
# Determine cell size based on the spacing of grid points
cell_size_x = (x.max() - x.min()) / (target_points_x - 1)
cell_size_y = (y.max() - y.min()) / (target_points_y - 1)

# Automatically determine the upper-left corner
upper_left_x = x.min()
upper_left_y = y.max()

# Create an affine transform for the grid (Y cell size is negative to go downward)
transform = rasterio.transform.from_origin(upper_left_x, upper_left_y, cell_size_x, -cell_size_y)
In [31]:
data_array = mag_grid
dMx, dMy, dMz = ev.compute_derivatives(data_array, transform, nodata_value, name, filename, mask, crs = None)
In [32]:
ev.plot_contours(x_grid, y_grid, dMx, label = 'nT/m', title = 'Horizontal Derivative in X-direction')
No description has been provided for this image
In [33]:
ev.plot_contours(x_grid, y_grid, dMy, label = 'nT/m', title = 'Horizontal Derivative in Y-direction')
No description has been provided for this image
In [34]:
ev.plot_contours(x_grid, y_grid, dMz, label = 'nT/m', title = 'Vertical Derivative')
No description has been provided for this image
In [35]:
max_depth = 1000
window_size=3
slide_distance=1
structural_index=1 
min_denominator=1e-8
In [36]:
distance_threshold = 7000
residual_threshold = 20000
# Apply Euler Deconvolution using a sliding window in X and Y directions
structural_index = 1  # Example structural index
start_time = time.time()
solutions = ev.euler_deconvolution_3d_window_xy_real(df, mag_grid, dMx, dMy, dMz, structural_index, window_size,
                                               distance_threshold, residual_threshold)
elapsed_time = time.time() - start_time
print(f"Time taken: {elapsed_time:.2f} seconds")

Processing Euler Deconvolution: 100%|██████████████████████████████████████| 39204/39204 [00:07<00:00, 4909.68window/s]

Time taken: 7.99 seconds
In [37]:
# Convert results to a DataFrame for easier visualization and further processing
solutions_df = pd.DataFrame(solutions)
In [38]:
solutions_df.dropna(inplace=True)
solutions_df
Out[38]:
x0 y0 z0 residuals
30 331750.0 6632600.0 210.067364 566.374462
94 335700.0 6632600.0 236.465625 1715.969104
95 336250.0 6633600.0 202.628918 750.396493
134 336950.0 6630000.0 195.879740 4.289047
137 337300.0 6629400.0 200.718215 2678.112174
... ... ... ... ...
30281 331250.0 6643000.0 138.670172 48.416218
30875 333350.0 6640400.0 171.739162 357.223035
32054 332850.0 6642800.0 100.257713 571.068698
35859 336200.0 6639000.0 191.712790 0.730324
36998 335050.0 6642400.0 106.255831 415.376270

10589 rows × 4 columns

In [39]:
sns.histplot(data=solutions_df, x='residuals', bins=100)
plt.xlabel('Residuals')
plt.ylabel('Count')
plt.title('Histogram of Residuals')
plt.show()
No description has been provided for this image
In [40]:
solutions_df = solutions_df[solutions_df['z0'] >=5]
In [41]:
solutions_df = solutions_df[solutions_df['residuals'] <=0.5]
In [42]:
sns.histplot(data=solutions_df, x='residuals', bins=100)
plt.xlabel('Residuals')
plt.ylabel('Count')
plt.title('Histogram of Residuals')
plt.show()
No description has been provided for this image
In [43]:
ev.plot_results(solutions_df, mag_grid, x_grid, y_grid)
No description has been provided for this image